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I. INTRODUCTION 

The theoretical description of the phase separation of 
fluid mixtures is a long-standing problem with impor- 
tance in many areas of soft matter physics such as, e.g., 
the stability of molecular and colloidal solutions 1 , the 
interactions between nanoparticles and macromolecules 
including novel phenomena such as Casimir forces 2 , as 
well as interfacial and confinement effects occurring in 
the presence of surfaces. Indeed, even "simple" mix- 
tures consisting of spherical particles can display non- 
trivial phase behaviour including triple points and criti- 
cal end points not present for one-component systems 3 ' 4 . 
Correspondingly, even more complex behaviour is ob- 
served for particles with internal degrees of freedom such 
as magnetic particles 5 ' 6 and mixtures involving shape- 
anisotropic particles such as colloidal rods 7 . 

Besides purely numerical approaches such as Monte- 
Carlo or Molecular (Brownian) dynamics computer sim- 
ulations, classical density functional theory (DFT) has 
proved to be a very accurate tool for describing both the 
homogeneous phase behaviour of mixtures and, at least 
for simple models, also the inhomogeneous structure oc- 
curring at interfaces. The key quantity in DFT is the 
one-particle density, which is obtained through minimiza- 
tion of a grand canonical free energy functional corre- 
sponding to the microscopic Hamiltonian of the system 8 . 
In addition to yielding the equilibrium phase diagram 
and microscopic fluid structure, DFT techniques have 
been successfully used to calculate nucleation barriers 
for state points in the mctastable region of the phase 
diagrams 9-12 . Based on the success of these approaches, 
it seems very tempting to use DFT techniques also to 
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tackle the non- equilibrium dynamics of the phase sepa- 
ration, including the growth of nuclei, the actual nucle- 
ation pathway and coarsening processes during spinodal 
decomposition. Indeed, one motivation for studying the 
phase separation dynamics accompanying demixing tran- 
sitions is their important role in the context of pattern 
formation and self-organization 13 ' 14 . Traditionally, phase 
separation dynamics is studied using mesoscopic models 
involving equations of motion for coarse-grained order 
parameters 15 . The advantage of addressing these topics 
using DFT is that the latter allows one to establish the 
link between the macroscopic behaviour of the system to 
the microscopic Hamiltonian, which is naturally incor- 
porated via the excess Hclmholtz free energy functional. 
In the last few years, the first steps in these directions 
have already been made on the basis of the so-called dy- 
namical density functional theory 16-19 (DDFT), which 
consists of a generalized continuity equation for the one- 
body density distribution of a many-particle systems of 
overdamped (Brownian) colloidal particles. Recent appli- 
cations of the DDFT to phase separation kinetics include 
spinodal decomposition in spherical fluids 18 and hetero- 
geneous nucleation at solid surfaces 20-22 . However, most 
of these studies have been devoted to simple fluids with 
no internal degrees of freedom. 

In the present work, we use both static and dynamic 
DFT to explore the phase separation of one of the sim- 
plest examples of a mixtures with internal degrees of 
freedom, that is, a binary fluid of spherical particles 
where one species carries a classical, (3D) Heisenberg 
spin. Heisenberg fluids 23-27 are basic models for contin- 
uum systems exhibiting ferromagnetic order, particularly 
for the description of ferromagnctism in undcrcooled liq- 
uid metal alloys. Mixtures of such systems and, in par- 
ticular, mixtures of magnetic and nonmagnetic particles 
are promising candidates for the controlled fabrication 
of patterns on the micron scale 5 . Moreover, an obvious 
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attractive feature of these systems is that the phase sepa- 
ration and thus, the occurrence of patterns can be tuned 
by external magnetic fields. 

From the theoretical side, the equilibrium properties 
of one-component Heisenberg fluids 23-27 , as well as other 
spin fluids with two-dimensional (XY) and Ising spins 
and mixtures thereof have been extensively studied 
by MC simulations, integral equation methods, and 
(mean field) density functional theories (see Ref. 28 
and references therein). However, this microscopic level 
of description for the dynamic behaviour is essentially 
unexplored. As a starting point to fill this gap we 
consider here a Heisenberg mixture in two spatial 
dimensions in the absence of an external field. The 
restriction to a 2D situation is actually close to many 
experiments (see, e.g., Ref. 5) and has the advantage 
that the calculated structure can be easily visualized. 
Within the vast parameter space characterizing our 
model, we focus on a situation where the system demixes 
into a non-magnetic and a ferromagnetic phase. For 
this situation, we first use conventional (static) DFT 
to calculate a complete phase diagram (involving a 
first-order transition and a tricritical point), as well 
as the inhomogeneous fluid density and magnetization 
profiles characterizing the liquid interface. Based on 
this information we then consider the phase separation 
dynamics, focussing on the nucleation of non-magnetic 
bubbles within the ferromagnetic liquid phase. For 
this problem, wc compare the results of three different 
approaches, namely classical nucleation theory (CNT), 
which is based on macroscopic concepts, an approach 
based on equilibrium DFT, and finally DDFT. Wc 
demonstrate that both DFT approaches yield consistent 
results for the nucleation barriers, but predict different 
pathways due to the fact that the DDFT conserves 
the densities (contrary to DFT). Moreover, the DFT 
nucleation barriers differ from the CNT predictions 
when the size of the critical nucleus becomes small. 
We also present evidence that the DDFT can describe 
the coarsening process during spinodal decomposition. 
The remainder of this paper is organized as followed. 
In Sec. II we formulate the model Hamiltonian for the 
binary system. The equilibrium theory of the demixing 
transition is presented in Sec. Ill, which includes a 
calculation of the phase diagram and the interfacial 
structure. In Sec. IV we turn to discuss the dynamics 
of the demixing transition. Finally, we summarize the 
results in Sec. V. 



II. MODEL 

The fluid system that we investigate is a binary mix- 
ture composed of two species. One species (A) is com- 
posed of spherical particles which interact via purely 
isotropic and repulsive forces. The other species (B) con- 



sists of magnetic particles. In addition to the repulsive 
interaction induced by the particle cores, these particles 
carry magnetic moments. The interaction part of the 
Hamiltonian may therefore be decomposed into a core 
part and a contribution from the spin-spin interaction 

where a, (5 = {A, B} and 

V a/3 {r l , Tj,8i, Sj) = V core (ri,r J -) + V r mag (r i , r j; s u Sj)5 a . B 5i3.B- 

(2) 

In our model, the particles are confined to a (two- 
dimensional) plane, so that the position of particle i, 
denoted = (xi,Z{), but the magnetic moment is repre- 
sented by a three-dimensional normalized classical spin 
Si whose orientation is described by the Euler angles 
lu = (0,(p). For the magnetic interaction, we choose the 
Heisenberg model, 

y mag (|r-r / |, W ,a; / )=J(|r-r / |) Sl .s 2 , (3) 

where J(|r — r'Q determines the range of the spin-spin 
interaction. We further assume that J(|r — r'Q can be 
described by Yukawa's potential, that is, 

(0, if |r - r'| < cr, 

J(|r-r'j)=<j T e -(|r-r'|/ g -p ( 4 ) 

I J |r-r'|/<7 ' ^ 

For interparticle distances |r — r'| < a the interaction be- 
tween two magnetic particles is assumed to be small as 
compared to the repulsion from the core potentials (sec 
below) and we therefore set the Yukawa potential in our 
model to zero in this region. The sign of the coupling 
constant J* = J/(k B T) in Eq. (4) (where k B is Boltz- 
mann's constant and T is the temperature) determines 
which type of spin ordering is preferred. As we show be- 
low in Sec. IIIB, the choice J* > yields a spontaneous 
ferromagnetic ordering of the magnetic component of the 
mixture at temperatures T below a Curie temperature 
Tc- On the other hand, the choice J* < favors anti- 
ferromagnetic ordering. To model the repulsion between 
the particles, we choose a Gaussian with height e and 
width a. The resulting "Gaussian core" model (GCM), 
first studied by Stillinger 29 , is given by 

W|r - r'|) = £cxp(-|r - r'| 2 /a 2 ). (5) 

The GCM is often used as an approximation for 
the effective interactions between the centre of mass 
of two "soft" particles, such as polymers and star- 
polymers 30 ' 31 or dendrimers 32-34 . The dimensionless 
quantity e* = e/{ksT) > determines the strength of 
the repulsion, and the range parameter a roughly corre- 
sponds to the radius of gyration of the 'particles'. For 
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the magnetic interaction, we choose a positive coupling 
constant J* > 0. Hence, in the ferromagnetic phase 
the magnetic contribution to the pair potential acts 
effectively as an attractive tail to the repulsive (soft) core. 



where A a denotes the thermal de Broglie wavelength of 
species a. The particle interactions enter into the excess 
part of the free energy functional, which can be written 

as 8 



III. EQUILIBRIUM THEORY OF THE DEMIXING 
TRANSITION 

A. The density functional 

The central quantity in density functional theory is the 
singlet (one body) density distribution p a (r, w). Follow- 
ing other studies of molecular magnetic fluids 23 we as- 
sume that the singlet density can be factorized into a 
translational (number density) part, p Q (r), and an orien- 
tational distribution function, h a (r, w), that is, 

p a (r,u) = p a (r)h a (r,uj). (6) 

The orientational distribution is normalized, so that, 



dujh a (r,Lu) = 1, 



(7) 



which yields Ha = 1/ (47r) for the particles without orien- 
tational degrees of freedom (species A). The equilibrium 
fluid singlet density distribution is that which minimizes 
the grand free energy functional 8 



E 



du \ dr 



V e xt(r,w) p a (r,u), 
(8) 



where T is the Helmholtz free energy functional and p a 
denotes the chemical potential for species a. Note also 
that Jdr = Jdx Jdz denotes a two-dimensional spatial in- 
tegral. For a given external potential T4 xt and interaction 
potentials, the functional in Eq. (8) has a minimum at the 
equilibrium density p° a {r,w). The functional Qy [{Pa}] 
is then identical to the grand canonical potential f2 (cf. 
Ref. 8). In the present study we set the external poten- 
tial T4 x t = 0. The Helmholtz free energy functional can 
be split up into two contributions: 



(9) 



where Jid[{Pa}] = E / dv I dujp a (r, cj)[ln(A 2 i o ct (r, w)) - 

a 

1] is the ideal gas contribution and .F ex [{/9 Q }] is the excess 
part. Using Eqs. (6) and (7), the ideal gas contribution 
becomes 



•Fid[{p Q }] =k B T drp A (r) \n(A 2 A p A (r)) - InAw - 1 



k B T J drp B {v) [ ln(A|p B (r)) - 



(10) 



+ k B T I drp B (r) / duh B (r, w) In h B (r, w)J , 

(11) 



^"ex[{Pi}] =2 E 



d\ dr dr' du doj' 



xp™(r,r>,u>,u/-,\)V a f > (\r-r'\,u>,u/). 



(12) 

Equation (12) is exact for systems with pair interactions. 

(2) 

The function p a p(r, r ,u, lo ; A) is the two-body density 
distribution function which is determined by the pair 
interactions y a ^(|r — r'\,u),w'), and A is a "charging" 
parameter 8 . In general, the function p^l is not known 
exactly. Here we employ a mean-field (MF) approxima- 
tion by setting /^(r, r', u, J; A) = p a (r, uj)p [3 {r', uj'), i.e. 
that the pair correlation function is set to one. Previous 
studies 31,35 have shown that this simple MF approxima- 
tion for the GCM yields reliable results for the fluid struc- 
ture and thermodynamics, particularly at higher densi- 
ties. Within the MF approximation, the excess part of 
the Helmholtz free energy may be written as the following 
sum 



•Fe: 



AA 



BB 



- T. 



AB 



BA 



(13) 



where the contribution due to the interactions between 
the non-magnetic particles is 



u r 



AA 



[{Pa}] = \ jdrJdr'p A {r)p A {r')V a 



(14) 



and for the magnetic particles 

^ B [{pb}] = \ JdwJckv'h B (r,uj)h B (r',uj') J dr J dr' 



x p B (r)p B (r') 



Korc(|r-r'|) + J(|r-r'|)s-s' . 

(15) 



The contributions T^ 3 and to the excess free en- 

ergy functional arc equal because of the symmetry of the 
pair potentials between species A and B (cf. Eq. (2)), 
that is, 



^ B [{ Pa }} = \ jdrJdr'pA(r)p B (r')V cole (\r-r'\) 



= [{ Pa }} 



(16) 



The equilibrium densities pa\r,u>) are found by mini- 
mizing the grand free energy functional given in Eq. (8): 



dTL[p a , h] 



dp a {r) 
dil[p a7 h] 
dh a (r, u) 



P< 0) (r) 



= o, 



0. 



,(0) 



(17) 



(r,u>) 
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Equation (17) yields an implicit equation for /ie(r,u;), 

tl B {r,UJ) — -y- — y -r y—r, (18) 

J c(cjexp(B(r) • s(w)) 
where the effective field 

B(r) = -[dr'[ du' p B {v')h B {v' , J) J(|r - r'|)s'. (19) 



Thus, the orientational distribution is determined solely 
by the scalar product between the spin and the effective 
field B(r). This is an exact result within the MF ap- 
proximation, which has been previously applied also to 
three-dimensional Heisenbcrg fluids 23 ' 26,27 as well as in 
other contexts such as in liquid crystal theory. 36 

In the low-temperature ferromagnetic state, the ori- 
entational order is uniaxial with respect to a director 
n. Thus, the angular distribution reduces to h B {r, w) = 
h B (r,u), where u = s ■ n = cos#. The integration over 
orientation in Eq. (19) then becomes 23,27 



JdLu'h B (r',u/)s' — > 2n J duh B (r, u)u = L(r), 



(20) 



where L(r) = coth_B(r) — 1/B(r) is the Langevin func- 
tion. The latter also defines the local magnetization 
m(r) = J duh B (r, u) cos 9 = L(r). Inserting Eq. (20) 
into Eq. (15), the Heisenberg contribution to the free en- 
ergy functional can be written as 

F™\p]=\ jdvjdv'p B {v)p B {v') 

x Kore(|r-r'|) + L(r)L(r')J(|r-r'|)". 

(21) 



B. Phase behavior 

In this section, we employ the density functional ap- 
proach introduced above to investigate the phase behav- 
ior of the bulk binary mixture. We restrict ourselves to 
fluid phases. The state of the system can be character- 
ized by the strength of the repulsion e* = e/(k B T), the 
magnetic coupling parameter 5 — J/e, the reduced total 
number density pa 2 in the system, and the concentration 
x of the magnetic component. The densities of the indi- 
vidual component can then be written as pa = (1 — x)p 
and p B = xp. Assuming that the system is homogeneous, 
the Helmholtz free energy per particle / = F/N follows 
from Eqs. (9)-(12) as 

f MF (p,x) = f id (p,x) + f™ F (p,x) 



Inx + (1 — x) ln(l — x) + x Jdu}h B (u})\n h B {uj) 

(1 - x) 2 V AA {0) + 2x(l - x)V AB (0) + x 2 V BB (0) 

(22) 



In Eq. (22), the coupling matrix elements V AA (0) and 
V AB (0) denote the integrated strength of the repulsive 
core potential, or equivalently the k — > limit of the 
Fourier transform of Eq. (5), that is, 



V AA (0)= drV coie (\v\) = V AB (0). 



(23) 



The element ^^^(O) involves, in addition, an integral 
over the magnetic interaction, i.e. 



V BB (0) 



drV c< 



M) 



+ 



drjdujj du'h B (u))h B (Lj')J(\r\)s ■ s' . 



(24) 



In the following, we investigate the possibility of fluid- 
fluid phase separation. The thermodynamic stability 
conditions for a binary mixture are given by 35,37 

(S).>* 

dx 2 



>0, 



a 2 / 

dv 2 



dx 2 



d 2 f 

dvdx 



>0, 



(25) 



where v = 1/p is the volume per particle. The first sta- 
bility condition expresses that the compressibility must 
be positive, the second ensures stability against sponta- 
neous demixing at constant volume, and the last inequal- 
ity is the condition for stability at constant pressure. As 
shown in previous studies 37 , it is more convenient to use 
these stability conditions in a constant-pressure ensem- 
ble. To this end, we perform a Legendre transform of the 
Helmholtz free energy per particle, yielding the Gibbs 
free energy per particle 



g{x,P) = f(x,v) 



df 
dv 



(26) 



For two phases I and II to coexist in equilibrium, the 
chemical potentials of each species a = A,B have to 
be equal, and the same holds for the pressure and the 
temperature. In other words, one has ^ a {p\,P B ) = 
MMyj), P 1 {A^b) = P ll (P l lP l B) and T\p\,p B ) = 
T ll {p\ l , p l g), where p^ 11 ^ arc the densities of the com- 
ponents in the two phases. These equilibrium conditions 
lead to a common-tangent construction on the Gibbs free 
energy, 



d_g_ 

dx 



dg 

dx 



g(xi,P) - g{x u ,P) 



xi - x u 



(27) 



where P is the (bulk) pressure at coexistence. The spin- 
odal is given by the inflection points of g, that is 



\dx 2 ) p 



0. 



(28) 
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FIG. 1: The phase diagram for a two-component GCM mixture, where one component is characterized by an 
additional Heisenberg interaction. The coupling parameters are a) e* = 5.0, 5 = 0.1 and b) e* = 5.0, 5 = 0.06. The 
yellow lines are isobars with line-to- line pressure difference of AP* = APa 2 / (fcgT) = 40. The black circles denote 
coexisting state points for P* — 40 (a) and P* = 80 (b), respectively. In a) the highest isobar indicated corresponds 
to P* = 160, and the (tri)critical point occurs at p c a 2 — 1.5, x c = 0.63. The corresponding data in b) are P* — 280, 

p c a 2 = 2.6, and x c = 0.61. 



Figure 1 shows the phase diagram for the (2D) bulk 
binary mixture whose free energy is given by Eq. (22). 
We consider a fixed repulsion strength e* = 5.0 and two 
different magnetic coupling parameters, S — J* /e* = 0.1 
and 5 = 0.06. For both parameter sets we find a first- 
order demixing phase transition appearing at densities 
above a critical density p c . Moreover, the demixing 
is coupled to a transition from a paramagnetic phase 
(in = 0), which is rich in A-particles, to a ferromagnetic 
phase (to > 0) rich in B-particles. The onset of magnetic 
order is determined by the Curie line (blue line in Fig. 
1), which is obtained by making a Taylor expansion of 
the Langevin function L [see Eq. (20) and below] around 
to = combined with the expression for the effective field 
B [see Eq. (19)] (note that we consider spatially homoge- 
neous systems here). The resulting critical concentration 
^Curic as function of the total density is given by 



•^-Curic (P> J*) 



2TrJ*pa 2 



(29) 



Inspecting the position of the Curie line in the phase di- 
agrams in Fig. 1, we see that the system is entirely dis- 
ordered, regardless of the concentration, for values of the 
total density pa 2 < 1 (1.5) for S = 0.1 (0.06). Increasing 
the density from these values towards the critical density, 
the transition from the paramagnetic into the ferromag- 
netic phase is of second order, until the Curie line meets 
with the demixing coexistence curve. This merging oc- 
curs directly at the demixing critical density (and crit- 
ical concentration), corresponding to a tricritical point. 
At densities p > p c , the Curie line then coincides with 
the low-concentration branch of the demixing spinodal. 
This reflects the fact that it is the spin-spin interaction 
[see Eq. (2)] which drives the phase separation. Indeed, 



as can be seen from Eq. (21), the spin-spin interaction 
reduces the free energy of the system whenever the mag- 
netization is non-zero. 

The demixing spinodal has been calculated using 
Eq. (28). Inside the spinodal the mixture is thermo- 
dynamically unstable, as indicated by the red areas in 
Fig. 1. The black areas in Fig. 1 indicate the metastable 
regions between the spinodal and the coexistence curve. 
The corresponding coexisting densities p l , p 11 and con- 
centrations x 1 , x 11 are calculated using Eq. (27). Recall 
that coexisting phases are at equal pressure; we display 
a number of isobars in the density-concentration plane, 
which arc indicated by the yellow lines in Fig. 1. For 
each value of the interaction parameter S we also dis- 
play, for one particular exemplary pressure value, a pair 
of coexisting state points (see circles). From this one 
may observe that the first-order phase transition is in- 
deed mainly a demixing transition in the sense that the 
total density change on crossing the transition is small. 
Finally, a comparison of Figs, la) and b) reveals that 
as S is reduced the critical point shifts to larger values 
of p. This is because an increase of the number density 
of the magnetic component supports the ferromagnetic 
phase transition [cf. Eq. (19)] . 



C. Interfacial structure 

In Sec. IIIB we showed that the magnetic mixture dis- 
plays a first-order demixing phase transition for a broad 
range of densities p and concentrations x. In this sec- 
tion, we focus on the structure of the fluid-fluid interface 
between the two demixed phases for densities above p c . 
The grand canonical free energy f2 for the non-uniform 
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binary mixture is given by Eq. (8). Setting the external 
potential V cx t — yields a well-defined (one-dimensional) 
interface between one region enriched with particles from 
the non-magnetic species and a second region enriched 
with magnetic particles. Setting the functional derivative 
of Eq. (8) to zero yields the Eulcr-Lagrange equations for 
the chemical potential of species A 



= 3 



p A = k B T\n p A (z) 



Pa(z') + pb(z') 



V' 



r'l), 



(30) 



and species B, respectively, 



p B = k B T{ln p B (z) +S(z)coth B(z) -1 




* 

p 


= 96.85 


* 

~ p 


= 71.47 


- p* 


= 62.39 


- p* 


= 53.96 



9 13 

Z to 



+ ln 



B(z') 



sinh [B{z')\ 



+ Jdv'[p A (z') + p B (z')\l3V coie (\r- r 
B(z) jdv'L B(z') PB (z')f3J(\r 



r'l), (31) 



where (3 = l/(k B T). Using the (bulk) coexisting densi- 
ties p A ulk = p\ , p B ulk = p l B found from the calculation of 
the binodal in section IIIB, 



the chemical potentials pa 
and p B can be eliminated from Eqs. (30) and (31). This 
leads to the equations for the one-body density profiles 
across the free interface. Specifically, we obtain for the 
non-magnetic component 




p « = 9685 
P* = 71.47 

P* = 62.39 

P* = 53.96 



z / o 



40 



PA(z) = p A 

B 

x exp 



£ (dv'[p^-p a {z'))pV cmc {\v-r'\) 

A J 



.a— A 

and for the magnetic component 

B 



(32) 



+ 



£jbulk£bulk _ B(z')L(z') 



Idv'(p^-P a {z'))^ 

*=A J 

+L{B h ^f |dr'(p^ lk - P 2(z'))pj(\r r'|) 

fi bulk B ( z ,j 

sinh(B bulk ) _ sinh[B(z'^ 
-l[b(z)] Jdr'L[B(z')] (^ ulk - p 2 (z'))f3J(\r r'\) 

(33) 

Equations (32) and (33) can be solved self-consistently. 
The results for the density profiles and the magnetization 
are shown in Fig. 2 for a fixed magnetic coupling param- 
eter 5 = 0.1 (see Fig. la) for the corresponding phase 
diagram). The decay of the density profiles into the two 
bulk phases in Fig. 2 is monotonic. In previous studies of 



r'l) 



FIG. 2: a) The density profiles of the two demixed 
phases and b) the magnetization profile (species B) at 
the interface for different values of P* . The coupling 
parameters are e* = 5.0 and 5 = 0.1. 



binary mixtures of soft particles 37,38 and also of colloid- 
polymer mixtures 39,40 it was found that non-monotonic 
oscillatory decay of the density profiles can occur for the 
free interface between coexisting state points that are suf- 
ficiently far removed in the phase diagram from the crit- 
ical point. We expect a similar scenario for the present 
system as that observed in the systems studied in Refs. 
37 and 38, i.e. we expect to observe oscillatory behav- 
ior in the density profiles also for the present system, 
but for higher values of the total bulk density, i.e. fur- 
ther away from the critical point. Note that the freezing 
transition 41 of the pure GCM fluid in three dimensions 
occurs only at much lower temperatures (i.e. much larger 
values of £*) than we consider here. 

Approaching the critical density from above, the in- 
terface softens, resulting in a reduced pressure parallel 
to the interface. Physically, this softening is reflected by 
a decreasing line tension (in 2D). To calculate the line 
tension we assume that the dividing surface is a line that 
is orientated along the x-direction. The pressure tensor 
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FIG. 3: The grand potential density as function of z 
with the interface being located at z = 25. 6(7 . The 
coupling parameters are e* = 5.0 and S = 0.1. 



FIG. 4: The interfacial (line) tension 7* = 7cr/(fc B T) 
between the demixed phases as function of the pressure 
difference relative to the critical point. The coupling 
parameters are 5 = 0.1, e* = 5.0. 



P is a (2 x 2)-matrix characterized by one component 
parallel to the interface (Pt) and one component normal 
to the interface (P/v), that is 



P T (z)e x e x + P N (z)e z e z 



(34) 



Here, e x and e z are normalized unit vectors in the in- 
direction and z-direction, respectively. In equilibrium, 
the normal pressure Pn{z) is constant and equals the 
bulk pressure P at coexistence. Furthermore, the tan- 
gential component Pt only depends on z. The interfa- 
cial (line) tension 7 is then defined as the excess force 
resulting from the dividing interface 42 ' 43 , that is, 



dz(P-Pr(z)). 



(35) 



42 



To actually calculate Pt(z), we use the relation 
Pt(z) = —w(z), where u(z) is defined as the grand 
canonical free energy per unit length calculated in the 
absence of an external potential. Numerical results for 
the function (oj(z)+P) [i.e., the integrand in Eq. (35)] at 
various total densities pa 2 are shown in Fig. 3. We only 
find non- vanishing values of the function (uj(z) + P) near 
the interface. This reflects the simple fact that the inter- 
facial tension stems from the density inhomogeneities at 
the interface (see Fig. 2). From Fig. 2 we see that the 
density profiles become smoother as the critical point is 
approached. This behavior is mirrored by the function 
(uj(z) + P) (see Fig. 3), resulting in a vanishing interfa- 
cial tension at the critical point. The behavior of 7 as a 
function of the pressure difference (P — P c ) is displayed 
in Fig. 4. 



IV. DEMIXING DYNAMICS 

We now turn to discuss the dynamics of phase sepa- 
ration in the present system. To this end, we employ 
the DDFT approach, 16-19 in which the time evolution 
of the one-particle densities are governed by a gener- 
alized continuity equation. The latter may be derived 
by integrating the Smoluchowski equation, that is, the 
Fokkcr-Planck equation for a system of (colloidal) parti- 
cles with overdamped stochastic equations of motion (i.e. 
the inertial terms in the microscopic Langevin equations 
of motion are neglected). The key approximation of the 
DDFT approach is that the non-equilibrium two-body 
density distribution functions at time t are set equal 
to those of an equilibrium system with the same one- 
body density profile 16-19 . As a consequence, the currents 
entering the DDFT equations are determined by (func- 
tional derivatives of) the equilibrium Helmholtz free en- 
ergy functional. As in most DDFT applications so far, we 
neglect here the effect of hydrodynamic (solvent-induced) 
interactions between the particles. 

For the present system, where one species (the B- 
particles) has internal degrees of freedom, one should 
employ the DDFT equations for anisotropic particles re- 
cently proposed in Rcf. 44, which shows that the equa- 
tions of motion for the position- and angle-dependent 
densities involve both the usual translational currents 
j a = —DaPaVSF/Spa, and also 'rotational current' 
terms resulting from application of the angular mo- 
mentum operator to the thermodynamic driving force 
5!F/5pB- However, the numerical solution of the re- 
sulting set of equations of motions for the demixing 
problems that we consider here (see below) involve si- 
multaneously determining not only the two-dimensional 
(isotropic) number density profiles of the A- and B- 
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FIG. 5: Density profiles of the isotropic particles (upper row), magnetic particles (middle row) and the 
magnetization (bottom row) as a function of the position. The time increases from the left to the right: t\ = 60tb, 
t 2 = 200t b , t 3 = 400r B and U = 800r B . The parameters are pa 2 = 3.2, x = 0.6, e* = 5.0 and 5 = 0.1. 



particles, respectively, but also of the orientational dis- 
tribution function /ig(r, The angle-dependence of 
the latter induces additional dimensions making the nu- 
merical calculations rather involved. 

In the present study, we somewhat simplify the prob- 
lem by assuming that the magnetic degrees of freedom 
are at each moment in time in equilibrium with the den- 
sity profiles. Physically, this assumption implies that 
the relaxation time of the magnetic moments is much 
shorter than that of the translational degrees of freedom. 
This assumption implies that the functional derivative 
8F '/6h(r,u),t) = at all times t, i.e., there is no driving 
torque. Under these conditions, the DDFT equations for 
the present system reduce to a coupled set of equations 
for the number density profiles of the A- and B- particles 



^dp a (r,t) 
dt 



= V 



Pa(r,*)V 



SJ r [p A (r,t),p B (r,t)} 



6p a (r,t) 



(36) 



combined with the self-consistency relation 
exp(B(r,i) • s(w,t)) 



J dui exp(B(r, t) ■ s(u), t)) ' 



(37) 



where the time-dependent effective field B(r, t) is given 
by Eq. (19) and where T is the MF Helmholtz free energy 
functional developed in Sec. Ill A. The mobility coeffi- 
cients r Q in Eq. (36) are related to the diffusion constants 
via T Q = D a j (fcflT), where a. = A or B. In what follows 
we assume that these are equal: Ta = Tb = T. 

When the external potential V cy ± = 0, the uniform 
density distributions, p a (r,t) = p„ ulk , always correspond 
to a stationary solution of the DDFT Eq. (36), since 
in this case the functional derivatives 5J-/Sp a are con- 
stants. However, if one considers applying small har- 
monic perturbations to the uniform densities p a (r,t) — 
Pa lllk + Sp a (r,t), where 5p a (r,t) ~ sin(k • r), with wave 
number |k| = k, and where 8p a a 2 <C 1, then one finds 
that inside the spinodal region (the red region in Fig. 1) 
density fluctuations with certain wavenumbers k grow 
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with time 18 ' 45 - i.e. within the spinodal the fluid is lin- 
early unstable. Outside the spinodal, the system is lin- 
early stable, i.e. the amplitude of any small amplitude 
density fluctuations decrease over time. In the region 
between the binodal and the spinodal (the black regions 
in Fig. 1), the fluid is linearly stable, but is not abso- 
lutely stable: if the amplitude of a given density pertur- 
bation is large enough, then the amplitude of this density 
fluctuation will grow over time. This is due to the non- 
linear terms in Eq. (36). For densities outside of the bin- 
odal, the uniform fluid is absolutely stable and all density 
modulations diminish in amplitude over time. The above 
description of the system is therefore qualitatively very 
similar to the results from Cahn-Hilliard theory 46 . 

Thus, there are two phase separation mechanisms: 
Firstly, spinodal demixing, which is triggered by the pres- 
ence of small amplitude density modulations, which in 
reality are always present due to thermal fluctuations. 
This is the dominant mechanism inside the spinodal. Sec- 
ondly, for state points in the region between the spinodal 
and the binodals, where the uniform fluid is linearly sta- 
ble, phase separation must proceed via the nucleation 
of density fluctuations with sufficiently large amplitude. 
We present results pertaining to these two mechanisms 
below. 



A. Spinodal demixing 

To study the phase separation dynamics in the spin- 
odal region, we set the time t = density profiles to be 
p a (r,t = 0) = p^ ulk + x( r ); where x( r ) is a small am- 
plitude random white noise field, which is equivalent to 
adding many harmonic density perturbations, with ran- 
domly chosen amplitude, phase and wavenumbers k. The 
density and magnetization profiles are then evolved for- 
ward in time using Eqs. (36) and (37). Note that we only 
add noise to the initial t = density profiles and do not 
add noise at any other subsequent time. This corresponds 
to taking a uniform system and then rapidly quenching 
it into the unstable region of the phase diagram, by de- 
creasing the temperature. Inside the spinodal, the den- 
sity perturbations with wave numbers < k < k c grow 
over time 18 ' 45 . The density modulations with wavenum- 
ber k rts fc*, where < fc* < k c , grows fastest in am- 
plitude over time, leading to density profiles having den- 
sity modulations with a typical length scale w 27r/fc >t , 
at short times after the quench. As is illustrated by 
the results displayed in Fig. 5, the domains of demixed 
fluid then coarsen over time. Since in the present system 
this demixing is between a magnetic and a non-magnetic 
phase, we also observe a similar pattern in the local mag- 
netization; see the bottom row in Fig. 5. The results in 
this figure are for a fluid with total density pa 2 = 3.2, 
concentration x = 0.6 and with e* = 5.0 and 5 = 0.1. 
The phase diagram for this system is displayed in Fig. 
1. For a given total density p within the spinodal region, 
on varying the concentration x one may observe bicon- 



tinuous labyrinthine patterns, such as those displayed in 
Fig. 5, or alternatively one observes phase separated mor- 
phologies consisting of 'islands' of the minority phase sur- 
rounded by a 'sea' of the majority phase. A more detailed 
discussion of spinodal phase separation and the resulting 
structures that we observe in the present system will be 
published elsewhere. 



B. Nucleation 

Before using DDFT to investigate the dynamics of nu- 
cleation, it is worth recalling the main results from clas- 
sical nucleation theory (CNT) and also approaches to 
nucleation using equilibrium DFT. 



1. Classical nucleation theory 

Nucleation is normally considered to be the relevant 
phase separation mechanism within the metastable re- 
gions of the phase diagram, i.e., within the regions be- 
tween the binodals and the spinodal. In these metastable 
regions, one may consider forming (circular, in 2D) clus- 
ters of the new (globally stable) phase with radius R, 
surrounded by the metastable bulk phase. One finds 
that the free energy as a function of R initially increases, 
reaches a maximum at R = R Cl the critical radius, and 
then decreases for R > R c . It is assumed that clusters 
of a given radius are randomly formed in the system by 
thermal fluctuations. Clusters with radius R < R c then 
typically shrink, since this lowers the system free energy. 
On the other hand, clusters with R > R c must grow 
without limit (in an infinite size system) , since doing this 
also reduces the system free energy, thereby initiating the 
transformation into the new (stable) phase. Thus, the 
excess free energy corresponding to the 'critical cluster', 
i.e. the cluster with radius i? c , is important, because this 
is the free energy barrier which must be surmounted for 
the phase separation to occur. 

CNT treats the nucleation process on the simplest pos- 
sible level. The key assumption of CNT is that any 
cluster (regardless of its actual size) can be regarded as 
a macroscopic object with a homogeneous density (and 
thus, pressure) inside and outside its surface. Moreover, 
interfacial curvature effects are typically neglected. As 
a consequence, the free energy related to creation of a 
nucleus can be written as a sum of two terms: a negative 
contribution stemming from the difference of the (bulk) 
pressures inside and outside the nucleus, and a positive 
contribution related to the increase of surface free energy. 
The latter is determined by the interfacial tension 7 of a 
planar interface (we use the interfacial tension calculated 
in Sec. IIIC above). Applying this concept to the two- 
dimensional system at hand, the grand potential for the 
creation of an ^4-rich phase nucleus in the surrounding 
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FIG. 6: The nucleation barrier Af2^ NT obtained from 

classical nucleation theory as a function of the 
concentration x for fixed total density per 2 = 3.2. The 
inset shows the excess grand potential of the nucleus as 
a function of the radius R for fixed concentration 
x = 0.96. The coupling parameters are e* = 5.0 and 
5 = 0.1. 



u sea" of S-rich phase is 



n 



CNT 



(R) = -P B A - ttR 2 \AP\ + 2irR~f, (38) 



where A is the total system area, ilo = —PbA is the 
grand potential of the uniform P-rich phase and AP 
is the difference between the bulk pressure Pb and the 
pressure of the A-rich phase that is being nucleated. 
The nucleation barrier is then given by the maximum 
of Afl CNT = fl GNT - fl which follows from Eq. (38) as 



Ail? NT = 



7T7 

IAPI 



The corresponding critical radius is 

7 



Rr 



IAPI 



(39) 



(40) 



For R > R c the grand potential decreases, indicating the 
absence of a sustaining force against the growth of a drop 
of the new phase. 

In the present calculations we focus on nucleation on 
the right-hand side of the phase diagram, i.e., the forma- 
tion of non-magnetic clusters predominantly containing 
A-particles, out of the ferromagnetic, P-dominated liq- 
uid. We consider the nucleus and its environment (i.e., 
the undersaturated magnetic liquid) at the same pair of 
chemical potentials fiA, P-B- Because we are working on 
the right side of the phase diagram, these chemical po- 
tentials are typically smaller than those corresponding 
to phase coexistence. The quantity AP appearing in 
Eqs. (38), (39) and (40) is then defined as the difference 
between the pressure corresponding to the actual (mag- 
netic) state on the right side, and the pressure of the 



FIG. 7: The grand potential Af2 of the nucleus 
obtained from the DFT approach as function of the 
excess number N cx for different concentrations. The 

parameters are pa 2 = 3.2, e* = 5.0 and S = 0.1. 



corresponding nonmagnetic state on the left hand side 
of the phase diagram. These states with equal chemical 
potentials (as well as the associated pressures) are found 
from the bulk free energy given in Eq. (22), which yields 
the desired quantities through the relations 



VA = f + P 
Pb = f + P 



dl 
dp 

df 



df 

dx 



P = P 2 



dpj x 



(41) 
(42) 
(43) 



In Fig. 6 we plot CNT results for the free energy bar- 
rier height Af2j? NT for a range of concentrations x of 
the magnetic species within the metastable regime. The 
chosen path corresponds to states at constant total den- 
sity per 2 = 3.2. The inset shows the grand potential [see 
Eq. (38)] as a function of R for one particular state point. 
As expected from the structure of Eq. (39), the nucleation 
barrier according to CNT becomes infinitely large at the 
coexistence line, where AP = 0. For the total density 
pa 2 = 3.2, the binodal is at the concentration x = 0.985. 
Decreasing then the concentration towards the spinodal 
value, the barrier height decreases. Note, however, that 
the barrier directly at the spinodal is not exactly zero, as 
one would expect at the limit of metastability. This de- 
ficiency is a well-known artefact of CNT. Indeed, given 
that CNT is a macroscopic theory, it is not surprising 
that its predictions become unreliable when the critical 
clusters become so small that they contain only a few 
particles. Under such conditions, a microscopic theory 
such as DFT is clearly more appropriate. 
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2. DFT approach to nucleation 

Since the important work of Oxtoby and Evans 9 there 
has been much work done using DFT to obtain a more 
reliable (microscopic) estimate for the free energy bar- 
rier Af2 c - see e.g. Refs. 10 ~ 12 ' 47,48 for examples of recent 
work. Here, we use the Oxtoby-Evans approach to study 
the nucleation of isotropic clusters of the phase rich in 
A-particles, and compare the resulting free energy bar- 
rier with the corresponding CNT results discussed in Sec. 
IV Bl. 

The key idea in all equilibrium DFT based approaches 
to nucleation is that the density profile characterizing the 
critical nucleus corresponds to a saddle point of the grand 
canonical free energy 9 . Assuming a symmetric droplet in 
the center of the system, the goal is thus to find the 
density profiles p a (r,cj) = p a (x,z 7 uj) for which 

6Q[{p a }] _ ST[{ Pa }} 



Sp a (r,uj) Sp a {r,oj] 



- Pa 



with the boundary conditions 



lim p a (x, z, u>) 

r— 7-00 



„bulk/ 



p~ a -'(x,z,uj) 

ra 



hulk h h a ulk (co), 



0, (44) 

(45) 
(46) 

where r = \J x 2 + z 2 . As demonstrated by Oxtoby and 
Evans 9 , these density profiles may be found by iterating 
the Euler-Lagrange equations, beginning with a simple 
approximation for the initial profiles. The latter are char- 
acterized by spherical symmetry and a sharp (step-wise) 
change of the number density at a radius R. The idea 
then is that if the guessed radius is too small (large) , the 
droplet will shrink (grow) during the iteration procedure 
until the profiles eventually approach the density values 
corresponding to the globally stable (unstable) phase. 
However, to identify the critical droplet one must iter- 
ate the Euler-Lagrange equations a limited number of 
times (for each initial guess R). This procedure allows 
for an estimate of the grand potential f2 as function of 
R. The critical droplet then follows as the position of the 
maximum. Inspired by these ideas we have performed ad- 
ditional calculations, not with the DDFT (which will be 
discussed below), but with a simple, relaxational algo- 
rithm which is equivalent to the Oxtoby-Evans method. 
Importantly, this algorithm keeps the chemical potentials 
fixed (just as in the original Oxtoby-Evans work 9 ), while 
the number densities themselves are not conserved. Us- 
ing this algorithm we investigated the evolution of several 
initial profiles of the form 




p a (x,z) 



. , 0, if r < R 
m(x,z) = < „ 

in , else 



(47) 



The initial values for the densities inside (p l a ) and out- 
side (p^ 1 ) the nucleus, are set to the bulk densities de- 
termined by the chosen value of the chemical potentials. 



The magnetization for these state points is obtained self- 
consistently from Eq. (19). For each initial guess, we 
performed about 300 iterations of the density profiles 
and the corresponding orientational profile, keeping the 
chemical potentials fixed. In this way we obtained an es- 
timate of the excess free energy AO for each given cluster 
size. As it turns out, the actual density profiles charac- 
terizing the clusters are rather smooth, such that the def- 
inition of a radius becomes ambiguous. As an alternative 
'reaction coordinate', we thus consider the quantity 



dxjdz( PA (x,z)-pY k ). (48) 



which counts the number of ^4-particles in the cluster 
(recall that we are considering the nucleation of isotropic 
A-dominated clusters). Numerical results for the func- 
tions Ail(N C x) at various chemical potentials within the 
metastable regime are plotted in Fig. 7. More precisely, 
the chemical potentials considered are associated to con- 
centrations (of B-particles) between the spinodal and the 
binodal along a path with fixed total number density 
pa 2 = 3.2 (CNT results along this path are presented in 
Fig. 6) . As may be seen from Fig. 7, all the curves reveal a 
clear maximum and thus, a clearly identifiable nucleation 
barrier, the height and position of which increase upon 
increasing the associated (bulk) concentration x. This 
finding is fully consistent with our expectation that the 
nucleation barrier and the size of the critical droplet are 
smallest close to the spinodal and then increase mono- 
tonically upon approaching coexistence. From the posi- 
tions of the maxima in Ail(N CK ) we also see that typ- 
ical critical droplets contain between a few ten and a 
few hundred of particles, consistent with results of other 
DFT studies 12 ' 48 . Our data for the height of the nucle- 
ation barrier as function of the concentration are summa- 
rized in Fig. 8, where we have included the corresponding 
macroscopic (CNT) results from Fig. 6. The main differ- 
ence between the two approaches is that the microscopic 
DFT calculation yields, contrary to CNT, a vanishing nu- 
cleation barrier at the spinodal, as one should expect on 
physical grounds. On the other hand, approaching the 
binodal the two curves merge, reflecting the increasingly 
macroscopic character of the critical cluster. 



3. DDFT approach to nucleation 

Whilst it is clear that the above and other equilib- 
rium DFT based approaches are able to calculate the 
density profiles corresponding to the critical droplet, it 
is not clear whether the other density profiles (i.e. those 
not corresponding to the critical droplet) have any phys- 
ical significance. Of particular interest are the density 
profiles corresponding to the most likely pathway (MLP) 
that go up to and then descend from the critical droplet 
state 12,47-49 . These MLP profiles should correspond to 
what one would observe experimentally for nucleation in 
the system. 
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FIG. 8: The nucleation barrier Af2 c obtained from 
classical nucleation theory (black curve) and from the 

DFT approach (red curve) as function of the 
concentration x. The dashed curve denotes the state 
point of coexistence. The parameters are pa 2 = 3.2, 
e* = 5.0 and S = 0.1. 




FIG. 9: The grand potential Af2 of the nucleus as 
function of the excess number N cx . The black curve is 
obtained from the DFT approach and the red curve 
corresponds to DDFT calculations for initial 
configurations with radius R — 2.05a (a), R = 2.5a (b), 
R = 2.75cr (c), R = 3a (d), R = 3.6a (e), R = 4a (f). 
The parameters are pa 2 = 3.2, x — 0.96, e* — 5.0 and 
(5 = 0.1. 



Recently, Lutsko 49 argued that one should be able to 
determine the MLP using DDFT by initiating the sys- 
tem at the critical droplet density profiles. Since these 
profiles correspond to a saddle point, and in reality one 
is never able to initiate the system exactly at the saddle 
point, the density profiles evolve under the DDFT away 
from this point. There are two directions the system 
may move: firstly, towards a state corresponding to the 
drop disappearing, with the excess density being spread 
uniformly throughout the system and secondly, the drop 
may grow. 

The performance of DDFT in this context is not yet 
fully understood (see Ref. 49), contrary to with the static 



DFT approaches mentioned above. One issue in this 
context is the fact that, within the DDFT, the num- 
ber densities are conserved quantities [as reflected by the 
appearance of a divergence in front of the free energy 
derivative in Eq. (36)]. This is in contrast to traditional 
(static) DFT approaches towards nucleation, where the 
fixed quantity is the chemical potential(s). In view of 
this subtle point, and given the rather plausible results 
from the DFT approach described so far, it is an impor- 
tant question whether a different algorithm, and particu- 
larly the conserved dynamics implied by the DDFT [see 
Eqs.(36)-(37)], yields consistent results. We recall that 
the DDFT is constructed such that the density evolves 
towards a profile which minimizes the free energy (one 
can prove that the free energy always decreases or re- 
mains constant under the time evolution of the DDFT, 
unless the system is externally driven 45 ). Therefore, one 
would expect that the critical profile found in the ap- 
proach discussed above, also 'behaves' as a saddle point 
within DDFT calculations. To confirm this, we have 
performed a number of DDFT calculations with initial 
density profiles stemming from the DFT calculations de- 
scribed above (after 300 iterations). Some results of these 
calculations are illustrated in Fig. 9. Note that for the 
DDFT results we have used a slightly different definition 
of N cyi . The DDFT calculations are performed on a fi- 
nite size square area of length L = 64a, with periodic 
boundary conditions. We start with the drop located at 
the center and we define 



N K 



dx / dz (pa{x, z) 



corner 

Pa 



(49) 



where p%™ T = pa (0,0) = p A (0,L) = p A (L,0) = 
Pa{L, L) is the value of the density at the corners of the 
(square, periodic) system. During the initial stages of the 
evolution p^ rnor = p A ulk , and so the value of N ex remains 
constant due to the conserved dynamics. However, at 
later times p™ rner changes, and so N ex changes with time. 
This is because either excess density from the center of 
the system diffuses out to the corner as the drop disap- 
pears, or because as the drop grows it removes particles 
of species A from the surrounding fluid and so p™ rnor de- 
creases. In Fig. 9 the curves labelled by a,b,c correspond 
to calculations where the radius characterizing the ini- 
tial profile is smaller than that characterizing the critical 
droplet predicted by the previous DFT approach. On the 
other hand, the curves labelled by d,e,f have been started 
from "supercritical" clusters. In all cases, the DDFT al- 
gorithm evolves in the direction predicted by the previous 
free-energy approach. That is, when starting from a sub- 
critical or supercritical profile, respectively, the droplet 
vanishes or grows without restriction. We also see from 
Fig. 9 that the actual values of Af2 "on the way" towards 
the final state are strongly different from those predicted 
by the previous approach. The interpretation of this is- 
sue clearly needs further investigation (see also the dis- 
cussion in Ref. 49). Furthermore, in the limit N cx — > 
the DDFT curves a,b,c approach different values, which 
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FIG. 10: Density profiles of the isotropic particles (species A) as function of the position. The upper row shows a 
supercritical growing nucleus (i? mltlal = 3.6<r) for subsequent times. The bottom row shows a time sequence for a 
subcritical nucleus (i? lmtlal = 2.5<r). The time increases from the left to the right. Upper row: t\ = Qtb, £2 = 200tb, 
£3 = 600tb and t 4 = 2000tb. Bottom row: t\ — Otb, ti — 40tb, £3 = 80tb and £4 = 200t£. The parameters are 

per 2 = 3.2, x = 0.96, e* = 5.0 and 6 = 0.1. 



depend on the initial profile. This is a consequence of 
the fact that the initial profiles a-c correspond to differ- 
ent space-averaged densities. This stems from the fact 
that the DDFT conserves the total densities, and so the 
final state emerging for subcritical clusters does not nec- 
essarily have chemical potentials equal to those chosen 
in the previous DFT calculations. We conclude that, 
at least, the DDFT predictions for the critical nucleus 
are consistent with the traditional DFT theory. Typical 
density profiles illustrating the nucleation dynamics ac- 
cording to DDFT on a microscopic (space-resolved) level 
are shown in Fig. 10. Besides the shrinking/growth pro- 
cess, we remark in particular the diffuse character of the 
density profiles, which directly reflects the difficulty in 
associating a fixed radius to the instantaneous droplets. 

V. CONCLUDING REMARKS 

In this paper we have investigated the demixing phase 
transition of a binary mixture of magnetic and non- 
magnetic soft-core particles from both a static and a dy- 
namic point of view. To this end we have employed clas- 
sical equilibrium DFT as well as the recently developed 
DDFT 16-19 , a time-dependent extension of DFT involv- 
ing a generalized continuity equation for the density. 

Our analysis of the equilibrium phase behavior of the 
system shows that the magnetic (Heisenberg) interac- 
tion within one species is capable of inducing macro- 
scopic phase separation. More precisely, for sufficiently 



large ferromagnetic coupling strength we find a combined 
phase transition where the system both demixes and de- 
velops global ferromagnetic ordering in one of the phases. 
Depending on the total density, the transition may ei- 
ther be second or first order (in terms of the magnetiza- 
tion and composition, respectively), with the two regimes 
being separated by a tricritical point. Thus, the gen- 
eral topology of the phase diagram (in the fluid regime), 
particularly the appearance of a Curie line and a tri- 
critical point, coincides with that of 3D Heisenberg and 
XY-fluids 28 ' 50 . Based on DFT we have also calculated 
the structure of the fluid-fluid interface in the first-order 
regime, i.e. we calculate the density and magnetization 
profiles, as well as the resulting interfacial (line) tension. 
As expected, the latter vanishes upon approaching the 
tricritical point from the high-density side. 

We note that the present DFT results are based on a 
mean-field (MF) approximation for the excess free energy 
contribution from both the soft-core repulsion (where MF 
theory has already proven to be very accurate 35 ) and 
from the Heisenberg interactions. Thus, the quantita- 
tive reliability of our results remain to be checked against 
those of more refined free energy functionals and/or com- 
puter simulations. However, judging from previous theo- 
retical studies of 3D Heisenberg fluids 23 we would expect 
the MF results to be quite good. 

The second part of the paper has been devoted to the 
dynamics of the (first-order) phase separation, concen- 
trating mainly on the nucleation of non-magnetic clus- 
ters within the metastable ferromagnetic phase. We 
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have analyzed and compared results from three dif- 
ferent approaches to nucleation. The simplest one is 
CNT, a macroscopic approach, involving the bulk pres- 
sure and line tension. The resulting nucleation barri- 
ers show familiar behavior (as compared to other model 
systems 12 ' 48 ) and predicts that the barrier height in- 
creases as one approaches to the binodal. In this regime 
CNT is expected to give a good approximation for AO c . 
However, CNT incorrectly predicts that the barrier is 
finite at the spinodal. As a second step we have inves- 
tigated the nucleation using the microscopic DFT ap- 
proach originally proposed by Oxtoby and Evans 9 . In 
this approach one calculates the free energy as function 
of an appropriate "reaction coordinate" , which we set to 
the excess number of particles in the nucleus. The re- 
sults for the nucleation barrier are very similar to those 
from CNT when the fluid state point is near to the bin- 
odal. However, in contrast to CNT, the DFT predicts 
that the nucleation barriers vanishes upon approaching 
the spinodal, as one should expect on physical grounds. 

Finally, we have touched the issue of the nucleation 
pathway by comparing the results of the microscopic 
DFT with those from DDFT. Since both of these ap- 
proaches are based on the same free energy functional, 
they are both probing the same underlying free energy 
landscape. The choice of reaction coordinate to some 
extent determines what regions of this landscape are ac- 
cessed and what path to and from the saddle point are 
predicted. In addition, the 'dynamics' inherent in each 
approach also determines this path. The DDFT has a 
real (physical) dynamics, which conserves the total num- 
ber of particles in the system, whilst in contrast the 
Oxtoby-Evans DFT approach, which has a finite number 
of iterations starting from an initial guess, leads to hav- 
ing an effective (unphysical) dynamics, which does not 
conserve the number of particles. Thus, the DDFT is es- 
sentially a canonical theory, whereas the DFT approach 
is more grand canonical in character. As we have shown, 
these subtleties are irrelevant for the actual height of the 
nucleation barrier, but yield marked differences when we 
consider the evolution of profiles away from that of the 
critical nucleus towards equilibrium. To our knowledge, 
these two DFT approaches have not been directly com- 
pared before. What is clear is that more investigations 
must be done to understand in detail the role of DDFT 
for the nucleation pathway, and more importantly which 
approach is the relevant one for determining the true 
MLP. 

There are several directions which we believe require 
further investigation. First, as a direct extension of the 
present work, it would be interesting to consider in more 
detail the spinodal decomposition occurring in the unsta- 
ble range of the phase diagram. Our results in the present 
study already indicate that the DDFT is capable of de- 
scribing the coarsening process characterizing this regime 
in a qualitative way (see Sec. IV A). However, this issue 
clearly calls for a more systematic investigation includ- 
ing a (linear) stability analysis 18 and also an investiga- 



tion of dynamic correlations such as the time-dependent 
structure factor. This quantity may be obtained using 
recent extensions of the DDFT to determine dynamical 
two-particle correlation functions 51 . 

Another (maybe more technical) open point concerns 
the effect of the actual formulation of the DDFT. In the 
present study we have restricted ourselves to a simplified 
version where the spin degrees of freedom relax instan- 
taneously; therefore, we were essentially dealing with the 
DDFT equations for isotropic particles, and the magne- 
tization entered only via the selfconsistency relation for 
the (local) magnetization. However, to further proceed 
one should also explore results from the full DDFT for 
anisotropic particles, which involves an equation for the 
rotational dynamics on top of that of the translational 
dynamics 44 . A particular interesting question in this con- 
text concerns the role of different time scales of rotational 
and translational motion for the demixing dynamics. 

Finally, it should be interesting to extend the present 
study towards 2D (magnetic) systems exposed to pat- 
terned surfaces and, optionally, additional external mag- 
netic fields. In fact, a number of recent experiments 52 ' 53 
have addressed the question of phase behavior and trans- 
port of magnetic particles on complex substrates 54 . All 
these question can be, in principle, tackled by the 
present (D)DFT approach. Work in these directions is 
in progress. 
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